{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QqIu96YO9lb8"
   },
   "source": [
    "# Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_bJFitQy9lcG"
   },
   "source": [
    "## An inferential problem: Wage Gap based on Sex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CnD-7-vI9lcH"
   },
   "source": [
    "In the previous lab, we already analyzed data from the March Supplement of the U.S. Current Population Survey (2015) and answered the question how to use job-relevant characteristics, such as education and experience, to best predict wages. Now, we focus on the following inference question:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BCOUK20N9lcI"
   },
   "source": [
    "**What is the difference in predicted wages between male and female workers with the same job-relevant characteristics?**\n",
    "\n",
    "Thus, we analyze if there is a difference in the payment of male and female workers (wage gap). The wage gap may partly reflect discrimination against female workers in the labor market or may partly reflect a selection effect, namely that female workers are relatively more likely to take on occupations that pay somewhat less (for example, school teaching).\n",
    "\n",
    "To investigate the wage gap, we consider the following log-linear regression model\n",
    "\n",
    "\\begin{align} \\log(Y) &= \\beta'X + \\epsilon\\\\ &= \\beta_1 D + \\beta_2' W + \\epsilon, \\end{align}\n",
    "where $D$ is the indicator of being female ($1$ if female and $0$ otherwise) and the $W$'s are controls explaining variation in wages. Considering transformed wages by the logarithm, we are analyzing the relative difference in the payment of male and female workers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Zjm2h8rt9lcJ"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import sklearn.linear_model as lm\n",
    "import statsmodels.formula.api as smf\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import Pipeline\n",
    "import warnings\n",
    "# ignore potential convergence warnings; for some small penalty levels, tried out, optimization might not converge\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Q2Z_SC7a9lcK"
   },
   "source": [
    "## Data Analysis\n",
    "\n",
    "We consider the same subsample of the U.S. Current Population Survey (2015) as in the previous lab. Let us load the data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LYwPvICz9lcK"
   },
   "outputs": [],
   "source": [
    "file = \"https://raw.githubusercontent.com/CausalAIBook/MetricsMLNotebooks/main/data/wage2015_subsample_inference.csv\"\n",
    "df = pd.read_csv(file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9UpPVZnh9lcL"
   },
   "source": [
    "To start our analysis, we compare the sample means for the different sexes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fw1c4JHm9lcL",
    "outputId": "1660091e-6bf6-4da2-ff55-09f806770df9"
   },
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gNhXyu279lcM",
    "outputId": "d7b624d7-695a-40ef-ab0d-eb46f88f1109"
   },
   "outputs": [],
   "source": [
    "table = pd.DataFrame()\n",
    "\n",
    "cols = [\"lwage\",\"sex\",\"shs\",\"hsg\",\"scl\",\"clg\",\"ad\",\"ne\",\"mw\",\"so\",\"we\",\"exp1\"]\n",
    "\n",
    "table['Variable'] = [\"Log Wage\",\"Sex\",\"Less then High School\",\"High School Graduate\",\"Some College\",\n",
    "                      \"Gollage Graduate\",\"Advanced Degree\", \"Northeast\",\"Midwest\",\"South\",\"West\",\"Experience\"]\n",
    "\n",
    "table['All'] = df[cols].mean().values\n",
    "\n",
    "table['Male'] = df[df['sex']==0][cols].mean().values\n",
    "table['Female'] = df[df['sex']==1][cols].mean().values\n",
    "\n",
    "\n",
    "# Show results\n",
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "o9K4z3T59lcN",
    "outputId": "45b832d8-a277-4d06-baba-441d4398f7ad"
   },
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HX98B_G19lcN",
    "outputId": "2cbb61a5-4644-4164-8870-5af315d8cd88"
   },
   "outputs": [],
   "source": [
    "# print to Latex\n",
    "print(table.style.to_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Q5wPwQk39lcO"
   },
   "source": [
    "In particular, the table above shows that the difference in average logwage between male and female workers is equal to  0,038."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RGPNaEUC9lcO",
    "outputId": "9c882b3c-7895-4afa-c39b-a3b1e31fce71"
   },
   "outputs": [],
   "source": [
    "# compute difference\n",
    "table.loc[0, 'Female'] - table.loc[0, 'Male'] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "h9Ax06WW9lcP"
   },
   "source": [
    "Thus, the unconditional wage gap is about  $3,8\\%$ for the group of never married workers (female workers get paid less on average in our sample). We also observe that never married female workers are relatively more educated than male workers and have lower working experience.\n",
    "\n",
    "This unconditional (predictive) effect of sex equals the coefficient  $\\beta$  in the univariate ols regression of  $Y$  on  $D$ :\n",
    "\n",
    "$$\n",
    "log(Y)=\\beta D+\\epsilon.\n",
    "$$\n",
    "\n",
    "We verify this by running an ols regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Kn6HiZnY9lcP",
    "outputId": "da7f7e03-d4f4-42b5-b628-7e69c7d0ac8b"
   },
   "outputs": [],
   "source": [
    "nocontrol_fit = smf.ols(\"lwage ~ sex\", data=df).fit()\n",
    "nocontrol_est = nocontrol_fit.params['sex']\n",
    "nocontrol_se = nocontrol_fit.HC0_se['sex']\n",
    "\n",
    "print(f\"The estimated sex coefficient is {nocontrol_est:.4f} \"\n",
    "      f\"and the corresponding robust standard error is {nocontrol_se:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7fOUkoYS9lcQ"
   },
   "source": [
    "Note that the standard error is computed to be robust to heteroskedasticity.\n",
    "\n",
    "Next, we run an ols regression of  $Y$  on  $(D,W)$  to control for the effect of covariates summarized in  $W$ :\n",
    "\n",
    "$$\n",
    "log(Y)=\\beta_1 D + \\beta_2' W + \\epsilon.\n",
    "$$\n",
    " \n",
    "Here, we are considering the flexible model from the previous lab. Hence,  $W$ controls for experience, education, region, and occupation and industry indicators plus transformations and two-way interactions.\n",
    "\n",
    "Let us run the ols regression with controls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DvaliSU29lcQ",
    "outputId": "e87a76e8-a3ea-4e06-8bb8-300bc3fb8259"
   },
   "outputs": [],
   "source": [
    "# Ols regression with controls\n",
    "\n",
    "flex = \"lwage ~ sex + (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)\"\n",
    "\n",
    "#   Note that ()*() operation in formula objects in statsmodels creates a formula of the sort:\n",
    "#  (exp1+exp2+exp3+exp4)+ (shs+hsg+scl+clg+occ2+ind2+mw+so+we) +  (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+occ2+ind2+mw+so+we)\n",
    "#  This is not intuitive at all, but that's what it does.\n",
    "\n",
    "control_fit = smf.ols(flex, data=df).fit()\n",
    "control_est = control_fit.params['sex']\n",
    "control_se = control_fit.HC0_se['sex']\n",
    "\n",
    "control_fit.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JFtJiI1K9lcR",
    "outputId": "ff1d5082-7c58-4091-a97e-62dd1cc268d0"
   },
   "outputs": [],
   "source": [
    "print(f\"The estimated sex coefficient is {control_est:.4f} \"\n",
    "      f\"and the corresponding robust standard error is {control_se:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zgH-fOfc9lcR"
   },
   "source": [
    "The estimated regression coefficient  $\\beta_1≈-0.06955$  measures how our linear prediction of wage changes if we set the sex variable  $D$  from 0 to 1, holding the controls  $W$  fixed. We can call this the predictive effect (PE), as it measures the impact of a variable on the prediction we make. Overall, we see that the unconditional wage gap of size  4 \\% for female workers increases to about  7 \\% after controlling for worker characteristics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l25upD8N9lcS"
   },
   "source": [
    "Next, we are using the Frisch-Waugh-Lovell theorem from the lecture partialling-out the linear effect of the controls via ols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dz9XPQTD9lcS",
    "outputId": "6aaa32d0-2430-4195-bb3f-d19062062045"
   },
   "outputs": [],
   "source": [
    "# Partialling-Out using ols\n",
    "\n",
    "# models\n",
    "flex_y = \"lwage ~  (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)\" # model for Y\n",
    "flex_d = \"sex ~ (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)\"   # model for D\n",
    "\n",
    "# partialling-out the linear effect of W from Y\n",
    "t_Y = smf.ols(flex_y, data=df).fit().resid\n",
    "# partialling-out the linear effect of W from D\n",
    "t_D = smf.ols(flex_d, data=df).fit().resid\n",
    "\n",
    "# regression of Y on D after partialling-out the effect of W\n",
    "partial_fit = sm.OLS(t_Y, t_D).fit()\n",
    "partial_est = partial_fit.params['x1']\n",
    "\n",
    "print(\"Coefficient for D via partialling-out \" + str(partial_est))\n",
    "\n",
    "# standard error\n",
    "partial_se = partial_fit.HC0_se['x1']\n",
    "\n",
    "# confidence interval\n",
    "print(\"95% CI: \" + str(partial_fit.conf_int().values[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"The estimated sex coefficient is {partial_est:.4f} \"\n",
    "      f\"and the corresponding robust standard error is {partial_se:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f_NwUdf09lcS"
   },
   "source": [
    "Again, the estimated coefficient measures the linear predictive effect (PE) of  $D$  on  $Y$  after taking out the linear effect of  $W$  on both of these variables. This coefficient equals the estimated coefficient from the ols regression with controls."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iskkYD3W9lcS"
   },
   "source": [
    "We know that the partialling-out approach works well when the dimension of $W$ is low in relation to the sample size  $n$ . When the dimension of  $W$  is relatively high, we need to use variable selection or penalization for regularization purposes.\n",
    "\n",
    "\n",
    "In the following, we illustrate the partialling-out approach using lasso instead of ols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "N_ZBt97j9lcT",
    "outputId": "9d8ee2aa-9078-42df-cec4-5b8358ddaf6f"
   },
   "outputs": [],
   "source": [
    "# Partialling-Out using lasso\n",
    "\n",
    "# models\n",
    "X_d =  smf.ols(flex_d, data=df).data.exog[:, 1:] # exclude intercept so that lasso does not penalize it\n",
    "y_d = smf.ols(flex_d, data=df).data.endog\n",
    "\n",
    "X_y =  smf.ols(flex_y, data=df).data.exog[:, 1:] # exclude intercept so that lasso does not penalize it\n",
    "y_y = smf.ols(flex_y, data=df).data.endog\n",
    "\n",
    "# partialling-out the linear effect of W from Y\n",
    "# train model using Lasso with cross validation and variable normalization\n",
    "lasso = Pipeline([('scale', StandardScaler()), # standardize the variables\n",
    "                  ('lasso', lm.LassoCV())])\n",
    "t_Y_lasso = y_y - lasso.fit(X_y, y_y).predict(X_y)\n",
    "\n",
    "# partialling-out the linear effect of W from D\n",
    "# train model using Lasso with cross validation and variable normalization\n",
    "lasso = Pipeline([('scale', StandardScaler()), # standardize the variables\n",
    "                  ('lasso', lm.LassoCV())])\n",
    "t_D_lasso = y_d - lasso.fit(X_d, y_d).predict(X_d)\n",
    "\n",
    "# regression of Y on D after partialling-out the effect of W\n",
    "partial_lasso_fit = sm.OLS(t_Y_lasso, t_D_lasso).fit()\n",
    "partial_lasso_est = partial_lasso_fit.params[0]\n",
    "\n",
    "print(\"Coefficient for D via partialling-out using lasso \" + str(partial_lasso_est))\n",
    "\n",
    "# standard error\n",
    "partial_lasso_se = partial_lasso_fit.HC0_se[0]\n",
    "\n",
    "# confidence interval\n",
    "print(\"95% CI: \" + str(partial_lasso_fit.conf_int()[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oRmJu6nm9lcT"
   },
   "source": [
    "Using lasso for partialling-out here provides similar results as using ols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "UpbveucD9lcU",
    "outputId": "84a087bf-75fd-441f-9640-66173a78e0d8"
   },
   "outputs": [],
   "source": [
    "table2 = pd.DataFrame()\n",
    "\n",
    "table2['Model'] = [\"Without controls\", \"full reg\", \"partial reg\", \"partial reg via lasso\"]\n",
    "\n",
    "table2['Estimate'] = [nocontrol_est, control_est, partial_est, partial_lasso_est]\n",
    "\n",
    "table2['Std. Error'] = [nocontrol_se, control_se, partial_se, partial_lasso_se]\n",
    "\n",
    "# Show results\n",
    "table2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZQnds_-i9lcU",
    "outputId": "2b73723e-faf6-4a2e-c0dd-4042e39d09a3"
   },
   "outputs": [],
   "source": [
    "# print to Latex\n",
    "print(table2.style.to_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jb17Q6AN9lcU"
   },
   "source": [
    "It it worth to notice that controlling for worker characteristics increases the wage gap from less that 4\\% to 7\\%. The controls we used in our analysis include 5 educational attainment indicators (less than high school graduates, high school graduates, some college, college graduate, and advanced degree), 4 region indicators (midwest, south, west, and northeast); a quartic term (first, second, third, and fourth power) in experience and 22 occupation and 23 industry indicators.\n",
    "\n",
    "Keep in mind that the predictive effect (PE) does not only measure discrimination (causal effect of being female), it also may reflect selection effects of unobserved differences in covariates between male and female workers in our sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iJWVmq--9lcV"
   },
   "source": [
    "Next we try \"extra\" flexible model, where we take interactions of all controls, giving us about 1000 controls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "G7L90Bpa9lcV",
    "outputId": "77acf08e-fc75-4bae-c36e-beeabfec589a"
   },
   "outputs": [],
   "source": [
    "# model\n",
    "extraflex = \"lwage ~ sex + (exp1+exp2+exp3+exp4+shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)**2\"\n",
    "\n",
    "control_fit_extra = smf.ols(extraflex, data=df).fit()\n",
    "control_extra_est = control_fit_extra.params['sex']\n",
    "\n",
    "n = df.shape[0]\n",
    "p = control_fit_extra.params.shape[0]\n",
    "\n",
    "control_extra_se = control_fit_extra.HC0_se[0]*np.sqrt(n/(n-p))\n",
    "# crude adjustment for the effect of dimensionality on OLS standard errors, motivated by Cattaneo, Jannson, and Newey (2018)\n",
    "# for really correct way of doing this, we need to implement Cattaneo, Jannson, and Newey (2018)'s procedure.\n",
    "\n",
    "print(\"Number of Extra-Flex Controls \" + str( p- 1))\n",
    "\n",
    "print(\"Coefficient for OLS with extra flex controls \" +  str(control_extra_est))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partialling-Out using lasso\n",
    "extraflex_y = \"lwage ~  (exp1+exp2+exp3+exp4+shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)**2\" # model for Y\n",
    "extraflex_d = \"sex ~ (exp1+exp2+exp3+exp4+shs+hsg+scl+clg+C(occ2)+C(ind2)+mw+so+we)**2\"   # model for D\n",
    "\n",
    "# models\n",
    "X_d =  smf.ols(extraflex_d, data=df).data.exog[:, 1:] # exclude intercept so that lasso does not penalize it\n",
    "y_d = smf.ols(extraflex_d, data=df).data.endog\n",
    "\n",
    "X_y =  smf.ols(extraflex_y, data=df).data.exog[:, 1:] # exclude intercept so that lasso does not penalize it\n",
    "y_y = smf.ols(extraflex_y, data=df).data.endog\n",
    "\n",
    "# partialling-out the linear effect of W from Y\n",
    "# train model using Lasso with cross validation and variable normalization\n",
    "lasso = Pipeline([('scale', StandardScaler()), # standardize the variables\n",
    "                  ('lasso', lm.LassoCV())])\n",
    "t_Y_lasso = y_y - lasso.fit(X_y, y_y).predict(X_y)\n",
    "\n",
    "# partialling-out the linear effect of W from D\n",
    "# train model using Lasso with cross validation and variable normalization\n",
    "lasso = Pipeline([('scale', StandardScaler()), # standardize the variables\n",
    "                  ('lasso', lm.LassoCV())])\n",
    "t_D_lasso = y_d - lasso.fit(X_d, y_d).predict(X_d)\n",
    "\n",
    "# regression of Y on D after partialling-out the effect of W\n",
    "partial_lasso_fit_extra = sm.OLS(t_Y_lasso, t_D_lasso).fit()\n",
    "partial_lasso_est_extra = partial_lasso_fit.params[0]\n",
    "\n",
    "print(\"Coefficient for D via partialling-out using lasso \" + str(partial_lasso_est_extra))\n",
    "\n",
    "# standard error\n",
    "partial_lasso_se_extra = partial_lasso_fit.HC0_se[0]\n",
    "\n",
    "# confidence interval\n",
    "print(\"2.5 - 97.5 % CI: \" + str(partial_lasso_fit_extra.conf_int()[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EyElA7c49lcY",
    "outputId": "33ef18a0-7ba3-4f09-e059-dc44fb0c0e47"
   },
   "outputs": [],
   "source": [
    "table3 = pd.DataFrame()\n",
    "\n",
    "table3['Model'] = [\"full reg\", \"partial reg via lasso\"]\n",
    "\n",
    "table3['Estimate'] = [control_extra_est, partial_lasso_est_extra]\n",
    "\n",
    "table3['Std. Error'] = [control_extra_se, partial_lasso_se_extra]\n",
    "\n",
    "# Show results\n",
    "table3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1NEU_9xj9lcY",
    "outputId": "9f55c7b2-2f59-42a8-8673-ac6bb83174c3"
   },
   "outputs": [],
   "source": [
    "# print to Latex\n",
    "print(table3.style.to_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cEg8OVKk9lcZ"
   },
   "source": [
    "In this case $p/n = 20\\%$, that is  $p/n$  is no longer small and we start seeing the differences between unregularized partialling out and regularized partialling out with lasso (double lasso). The results based on double lasso have rigorous guarantees in this non-small $p/n$ regime under approximate sparsity. The results based on OLS still have guarantees in $p/n< 1$ regime under assumptions laid out in Cattaneo, Newey, and Jansson (2018), without approximate sparsity, although other regularity conditions are needed."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "PM1_inference.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
